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METHOD OF ESTIMATING THE ILLUMINATION FOLD IN THE 

MIGRATED DOMAIN 

The field of the invention is seismic prospecting. 
The invention relates to seismic imagery and more 
particularly to estimating the illumination fold in the 
migrated 3D domain necessary to obtain a high resolution 
5 image of the underground structure at true amplitude. 

Seismic prospecting generally consists of emitting 
seismic waves into the sub-soil using one or several 
seismic sources, and recording seismic data corresponding 
to seismic waves reflected on geological interfaces in 
10 the sub-soil (also called reflectors) on the ground 
surface, . as a function of time, using receivers (also 
called geophones or hydrophones depending on whether the 
survey is _ being made on land or at sea) and then 
processing these data to extract useful information about 
15 the geology of the sub-soil. 

Data processing (or seismic traces) recorded by 
receivers are useful particularly for designing images 
with information about the geology of the sub-soil. 

The objective of high resolution seismic imagery 
20 with true amplitude is to provide an image of the sub- 
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soil that is not disturbed, particularly by wave field 
propagation effects in the sub-soil and by acquisition 
effects (we then talk about acquisition patterns) . 

Acoustic reflectivity of the sub-soil is represented 
5 as an image using a so-called migration method to restore 
correct forms of geological interfaces. To achieve this, 
seismic data are converted from the domain in which they 
are acquired to the migrated domain in which geological 
interfaces are represented in their genuine locations. 
10 In the migrated domain, the vertical scale may be 

time (we then talk about time migration) or depth (we 
then' talk about depth migration) . 

Migration is usually done before summation (we then 
talk about pre-stack migration) , with direct mapping of 
15 the seismic signal from seismic traces towards the 
migrated domain . 

The precise identification of the sub-soil geology 
at a particular point is directly related to the 
capability of a seismic emission to illuminate this 
20 point. If a point in the sub-soil is highly illuminated, 
in all likelihood the seismic survey can evaluate the 
geology of the sub-soil at this point. On the other hand, 
if a point of the sub-soil is not well illuminated, it is 
quite possible that the seismic survey will provide an 
25 imprecise or even incorrect evaluation of the sub-soil 
geology at this point. 

In order to obtain a precise representation of the 
sub-soil, the true amplitudes of the reflection of. the 
seismic wave on a reflector must be represented 
30 perfectly. 

Therefore, it is necessary to know precisely if a 
change in the amplitude at a specific point in the sub- 
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soil is the result of a change in the geology of the sub- 
soil at this point or if it is the result of a. particular 
illumination at this point. 

But reflectors are usually illuminated non- 
uniformly. In particular, the illumination depends on: 

- the configuration of seismic data acquisition; 

- lateral and vertical variations of the propagation 
velocity of seismic waves in the sub-soil. 

In order to obtain a high resolution image of the 
reflectivity of the sub-soil, it is necessary to make an 
amplitude correction to take account of non-uniformity of 
illumination of reflectors when migrating 'seismic 
reflections . 

As has already been described above, this correction 
(done using a weighting factor) is necessary because: 

- the reflectivity of the sub-soil . cannot be 
precisely imaged without correcting amplitude distortions 
related to propagation of the wave field in the sub-soil 
(spherical divergence, absorption, directivity, etc.); 

- to improve the resolution of the reflectivity 
image, it is important to take account of the variability 
of the migration illumination fold for each inclination 
of a reflector used during the migration. 

For a given image point in the migrated domain, the 
illumination fold corresponds to the number of source- 
receiver pairs for which contributions for building up 
the amplitude of the reflectivity image interfere 
constructively (stationary phase condition) . 

The document "Bleistein, N. 1987, On the Imaging of 
Reflectors in the Earth: Geophysics Soc. of. Expl . 
Geophysics, 931-942" proposes a method used to obtain an 
amplitude weighting factor for a regular acquisition 
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configuration, capable of making a high resolution 
representation of the reflectivity of the sub-soil in 
true amplitude. 

This is an "inverse" method composed of a modelling 
5 step (modelling of seismic traces from a reflectivity 
model) followed by a migration step (from modelled 
seismic traces to imaged reflectivity) , the weighting 
factor then being chosen so as to obtain perfect 
agreement between the imaged reflectivity and the 
10 reflectivity model. 

By following this "inverse" method proposed by 
Bleistein and as shown in Figure 1, it can be seen that: 

- for a given acquisition configuration 
(distribution of (source S, receiver R) pairs), 

15 - for a given image point x in the migrated domain, 

- for a given reflector • inclination (indicated by a 
dip vector p ) , 

the illumination fold \(x,p) is determined by the 
number of source-receiver pairs such that their migration 
20 . contributions constructively interfere with the image of 
the dip reflector p at image point x . 

Expressions developed by Bleistein for the 
illumination fold \{x 9 p) include the effect of the entire 
geometric configuration of the acquisition in a single 
25 step. 

Consequently, the method proposed by Bleistein can 
only be used for a regular and parametric distribution of 
all source-receiver pairs. 

Consequently, this restrictive application condition 
30 limits its use to a small number of acquisition 
configurations. 
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Therefore for acquisition configurations with 
irregular geometries, it is necessary to define 
individual illumination folds | ( x,p,sj ) , for each (source 
s, receiver r) pair involved in the' migration. 
Obviously, the said individual ratios contribute to the 
total illumination fold |(3t,^) according to I{x,p) 
(x,p,s,r ) ) . 

s.r 

However, the study carried out by Bleistein is based 
on a number of developments that were intended to 
estimate the illumination fold in the migrated domain for 
any geometric acquisition configuration, and for any 
distribution of the source-receiver pairs. 

For years, various approximations have been proposed 
for the seismic wave that is, the key element for 
estimating the stationary phase condition (constructive 
interference condition) . 

The first developments consisted of modelling 
specular reflection points along a horizon (2D map) 
selected in the migrated domain (3D) by ray tracing for 
each of the (source s, receiver r ). pairs . 

This reflection modelling step determines the 
following quantities: 

- location of reflection points x r {s,r) along the 
selected horizon; 

- the specular reflection travel, time t r ( x r ; 5" , r ) ; 
The diffraction travel time t d along the selected 

horizon can then be estimated by a development of the 
reflection travel time t r , around the specular reflection.' 
point x r (s,r) in the Taylor series (determination of 
first and second derivatives V x t r (x r ) and A xx t r (x r )). 
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Figure 2 illustrates this procedure during which the 
diffraction .travel time t d ( x ) at image point 5c is 
estimated using the reflection travel time t r (x r ) and 
more particularly by its development around the specular 
ref lection point x r in the Taylor series. 

Summation of all wavelet mappings obtained for each 
source-receiver pair then provides a means of estimating 
the illumination fold in the migrated domain along the 
selected horizon. 

By modelling the seismic wavelet according to a 
pulse peak, with infinite frequency band, document "Cain, 
G., Cambois, G., Gehin, M and Hall, R., 1998 Reducing 
risk in seismic acquisition and interpretation of complex 
targets using a Gocad-based 3-D modelling tool, 68 th Ann. 
Internat. Mtg, Soc . of Expl . Geophys . , 2072-2075" has 
shown that an impact map can be obtained, and in this 
case only the specular contribution on a selected horizon 
is actually taken into account. 

Considering a realistic approximation of a seismic 
wavelet with limited band, the following documents: 

- "Schneider, W. , Jr and Winbow, G. , 1999, Efficient 
and accurate modelling of 3D seismic illumination, 69 th 
Ann. Internat. Mtg, Soc. Of Expl. Geophys, Expanded 
Abstract, 633-636"; 

- "Laurain, R. and Vinje, V. 2001, Pre-stack 
migration and illumination maps, 71th Ann. Internat., 
Mtg: Soc. Of Expl. Geophys, Expanded Abstract, 92 9-932"; 
have shown that it is possible to improve estimating the 
illumination fold along selected horizons. 

These approaches, in particular based on scanning of 
selected horizons, have been found to be reliable and 



WO 2005/008292 



7 



PCT/IB2004/002618 



precise. However, their efficiency is limited 
particularly due to the. fact that: 

- compensation for the illumination fold is limited 
to a finite number of selected horizons and therefore, 
unlike the method developed by Bleistein, this method is 
not completely 3D; 

- for each (source s , receiver r) pair involved in 
the migration, a dedicated modelling step must be carried 
out after the migration. However, this modelling step to 
determine the (x r9 V x t r9 A XtX t r ) quantities mentioned 
previously by ray tracing along each of the selected 
horizons is particularly expensive in terms of time. 

In order to overcome these limitations, and in order 
to approximate a 3D approach, a two-step solution has 
been proposed to estimate, and then to compensate for the 
illumination ratio, and particularly in documents: 

- "Albertin, Bloor, R. , Beasley, C, Chang, W., 
Jaramillo, H. , Mobley, E. and Yingst, D., 1999, Aspects 
of true amplitude migration, 69 th Ann . Internat . Mtg: 
Soc. Of Expl. Geophys, Expanded Abstract, 1247-1250"; 

- "Audebert, F. , Froidevaux, P. Huard, I., 
Nicoletis, L. and Svay-Lucas, J., 2000,- A multi angle 
toolbox for restored amplitude images and AVA gathers, 
7 0th Ann. Internat, Mtg: Soc. Of Expl. Geophys, Expanded 
Abstract, 1020-1023". 

The first step consists of making an analysis of the 
multi-dip illumination fold . using predetermined tables 
for Kirchhoff migration for diffraction travel times 
t d (x;s,r), in the migrated domain. 

In the second step,' the inverse values of the 
estimated -illumination folds are used to make the 
weighted summation of seismic data. 
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For each image point x and for an entire reflector 
dip range p , it has been proposed to increment an 
illumination occurrence counter I(x 9 p) every time that a 
migration contribution reaches the image point x with 
5 the dip component p . • 

In a manner similar to that proposed by Cain et . al 
(see above), the criterion used to increase the 
illumination occurrence counter takes account of a pulse 
model of the seismic wavelet, with infinite frequency 
10 band. 

This high frequency approximation is such that for a 
given dip component, only specular migration 
contributions corresponding to the said dip component are 
used for estimating the illumination. 
15 And there are a number of limitations to this high 

frequency approximation: 

- the estimated illumination is biased since the 
limited band nature of the seismic signal is ignored; 

- in order to avoid division by zeros, the 
20 illumination folds are stored by dip classes. However, 

the size of each class is a critical element. If the 
class size is too small, risks of numerical instability 
are very high; if the class size is too large, the 
precision of the estimated illumination coverage .is low. 

25 The invention is designed to provide a method for 

estimating the seismic illumination fold in the migrated 
3D domain that is not subject to these limitations. 

More precisely, the purpose of the invention is to 
propose a method that takes account of the limited band 

30 nature of the seismic signal in order to estimate more 
realistic illumination folds, and which at the same time 
solves the problem of the class size. 
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Consequently, according to a first aspect of the 
invention, it divulges a method for estimating the 
seismic illumination fold \{x,p) in the migrated 3D 
domain at an image point x for a dip of vector p 

5 characterised in that it estimates the illumination fold 
I(x 9 p 9 s,r) for each (source s, receiver r) pair in the 
seismic survey, by applying the following steps: 

- determination of the reflection travel time 
t r (x r (p);sr) from the source s to the specular 

10 reflection point x r on the plane reflector passing 
through the image point x and perpendicular to the dip 
vector p and then return to the reflector re- 
starting from the diffraction travel time t d (x;I,r) 
. from the source s to the said image point x and then 
15 return to the reflector r ; 

- incrementing the said illumination fold I (x, p, J, r) 
related to the said (source J, receiver r) pair as a 
function of the difference between the diffraction travel 
.time t 6 (x;s,r) and the reflection travel time 

20 t r (i r (p);fr). 

Preferred but non-limitative aspects of the method 
according to the first aspect of the invention are as 
follows : 

- the method also comprises a step to summate each 
25 illumination fold I (x, p, J S F) related to a (source s , 

receiver r ) pair so as to determine the total 
illumination fold I(x,p) = ^I(x 9 p 9 s 9 r) , 

s,p 

- during the increment step, the illumination fold 
I(x 9 p,s,r) may be' incremented using an increment function 
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i(t d ,t r ; s,r) using the equation I{x,p,s,r) = I(x,p,s,r) + 
i(t d ,t r ; s,r), 

- the increment function i may be a function of the 
seismic wavelet s(t), 

5 - the increment function i may be expressed as a 

function of the derivative of the seismic wavelet with 

respect to timer 

- it is envisaged that an a priori correction 
w(jc,j,r) of the illumination fold could be taken into 

0 account, which usually has to be considered during 
migration, 

- the determination step includes the second order 
Taylor series development of the diffraction travel time 
t d (x;s,r) around the image point x , 

5 - the specular reflection point x r {p ) is determined 

along the length of the said reflector such that the 
diffraction travel time at the said specular reflection 
point x r {p) is stationary, 

- the determination step uses isochronic migration 
20 maps t a (x;s,r) specified for each (source s- r receiver r ) 

pair involved in the migration and each image point x. in 
the migrated 3D domain, 

- the seismic illumination fold I(x,p) in the 
migrated 3D domain is estimated during the Kirchhoff 

2 5 summation migration of seismic data recorded during the 
3D seismic prospecting. 

According to another aspect, the invention proposes 
a method for correction of seismic data amplitudes 
recorded during a 3D seismic survey .in order to 

30 compensate for the effect of non-uniform illumination of 
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sub-soil reflectors, characterised in that if comprises 

steps consisting of: 

- estimating the illumination fold I(x,p) using the 

method according to the first aspect of the invention, 

- using the inverse I" 1 ^) of ' the said ratio as a 
weighting factor to be applied to each of the said 
seismic data amplitudes. 

According to another aspect, the invention proposes 
a method for selection of an acquisition geometry among a 
plurality of acquisition geometries as a function of the 
target of 3D seismic prospecting, characterised in that 
it consists of the following steps: 

- determining the illumination fold I(x,p) by the 
method according to the first aspect of the invention, 
for each of the acquisition geometries considered, 

- selecting the acquisition geometry providing the 
' optimum illumination fold as a function of the target. 

Other aspects, purposes and advantages of the 
invention will become clearer after reading the following 
detailed description of a preferred embodiment of it, 
given as a non-limitative example with reference to the 
attached drawings on which: 

- Figure 1, already commented, illustrates the 
constructive or non-constructive contribution as a 
function of s (t d -t r ) , of the seismic trace migrated to 
the image of the reflector with a dip along x ; 

- Figure 2, already commented, illustrates the 
estimate made according to the state of the art for the 
diffraction travel time at the image point t d ( * ) starting 
from a Taylor series development of the reflection travel 
time at the specular reflection point t r U r ); 
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- Figure 3 illustrates the method according to the 
invention according to which the specular reflection 
point x is estimated, together with the reflection 
travel time t r (*r) starting from a Taylor series 
development of the diffraction travel time at the image 
point t d ( x ) . 

- Figure 4 is a flowchart diagrammatically showing 
the main steps in the method according to the invention. 

The invention proposes a method of estimating the 
illumination fold in the migrated domain I(x,p;s,r) 
for any migrated position x, regardless of its dip 
component p and regardless of the distribution of 
seismic traces, using only the most basic quantities 
involved in migration by a Kirchhoff summation, namely: 

- the coordinates (s,r) of source-receiver pairs, 
where s represents the position of a source S and r 
represents the position- of a receiver R; * 

- migration isochronic maps td(x;?,r) for each of 
the (source s , receiver r) pairs involved in the 
migration. 

The invention is equally applicable to imagery in 
time and imagery in depth. 

The invention divulges a method of the type 
including a modelling step followed by a migration step. 

Concerning the modelling step of the reflection 
point x and the reflection travel time t E (x r ), the 
method according to the invention applies the stationary 
phase condition to a Taylor series development of 
diffraction travel times . t d ( x ) . 

The said development is a second order development. 

Concerning the migration step, the method according 
to the invention uses isochronic migration maps available 
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for the migration conventionally done by a Kirchhoff 
summation . 

For each of the (source s , receiver r) pairs, the 
illumination fold I ( x , p ; J , f ) is updated with reference 
to Figure 3 using the following equation: 

I(x,p;s,r) = Hx,p;s,r) + I(td,tr; s,r) Equation (1) 
where 

- x is the image point in the migrated domain x = 
(xi, x 2 , x 3 ) T ; 

- p is the dip vector representing the selected dip 

component, p = (Pi, P2» P3) ' 

- s defines the coordinates of the source, s = (s lf 

s 2 , s 3 ) T ; 

- f defines the coordinates of the receiver, r = 
(ri, r 2 , r 3 ) T ; 

- td(3c; s , r) is the diffraction travel time from 
the source J, as far as the. migrated position x, and 
return as far as the receiver f (path indicated by t d on 
Figure 3 and shown in solid lines); 

- t T (Xr(p); s, r) is- the reflection travel time 
from the source s as far as point x of specular- 
reflection on the plane reflector passing through image 
point x and perpendicular to the dip vector p, and 
return as far as the receiver r (path indicated by t r on 
Figure 3 and shown in dashed lines); 

- I(td,tr; s , r) is an increment function for which 
the expression will be described in detail below. 

As already mentioned above, the method according to 
the invention advantageously uses isochronic migration 
maps t d (3c; s , r) involved in the migration by Kirchhoff 
summation, for each of the (source s , receiver r) pairs 
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involved in the migration. These maps are "input data" to 
the method according to the invention. 

It also uses first and second derivatives 
V y t d (5c; s,r ), A x , x t d ( 3c ; s , r ) of migration isochrones when 
they are available. If not, the said derivatives are not 
input data, but may be calculated as part of the method 
according to the invention. 

Figure 4 is a flow chart showing the different steps 
in the method according to a preferred embodiment of the 
invention. 

The said method, consists of carrying out a step 1 
used for each image point x in the ' migrated 3D domain, 
for each of the (source s, receiver r) pairs involved in 
the migration. 

The said step 1 includes two operations 2, 3, for 
each dip component p selected for the illumination 
calculation I ( x , p ; s , r ) . 

The first operation 2 consists of determining the 
reflection travel time t r (x r (]5); s , r) to the specular 
reflection point xAp) on the plane reflector passing 
through the image point x and perpendicular to the dip 
vector p . 

Operation 2 is determined from the diffraction 
travel time td ( 5c ; s, r) known from the isochronic 
migration maps. 

The second operation 3 consists of incrementing the 
illumination fold in- the migrated domain l(x , p ; s , r ) of 
an increment function I (td, tr; $, r), the sard 
increment function taking account of: 

- the difference {t d (x;s,r) - t E ( x r ( p ) ; s , r ) } 
between the travel times mentioned above-; 
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- the seismic wavelet s(t); 

- and possibly the migration amplitude term 
w ( x ; s r ) . 

At the end of operations 2 and 3, the illumination 
5 amplitude ratio I{x,p;s,r) is known for a given image 
point x and a dip component p . 

By repeating the operations mentioned above for each 
of the (source s , ' receiver r) pairs, the final step is 
to determine the illumination ratio in the migrated 
10 domain I ( x , p ) . 

This is done by adding the relative illumination to 
each pair (source s, receiver r) according to I ( x , p ) = 

(x, p ; s r f ) . 

The reflection point x r and the reflection travel 
15 time t r are determined (during operation 2 mentioned 
above) starting from a stationary phase analysis, by a 
second order Taylor series development of the diffraction 
travel time t d around the image point x . 

The following expressions for x r and t r (equations 
20 (3) and (4) ) are obtained from this analysis .(see 
equation (2) below): 

t d (X;S , r )=t r (<x;J, r ) + ( V >: t r ) r ( x-x r ) + ~ ( x -x r ) T (V :li:: t T ) (x-5 r )+... 

Equation (2 ) 

25 x*{p) = Jt-M.F" 1 . b Equation (3) 

t r ( x r (P ) ; s , r ) = t d {x ; s , r ) ~~ . b r . F" 1 . b 

Equation (4) 

where, in equations (3) and (4). 
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- M is a (3x2) matrix described by two vectors 
extending along the length of the reflection plane. 
Consequently, these two vectors are perpendicular to the 
dip vector p . 

5 - b is a (2 x 1) vector of the second order 

derivatives of the diffraction travel time along the 
reflection plane, 

b = M T . (V x t d ) Equation (5) 

-F is a (2x2) matrix of first order derivatives 
10 of the diffraction travel time along the reflection 
plane, 

F = M T . (A XiX t d ) . M Equation (6) 

The following describes how equations (3) and (4) 
are obtained. 

15 Note that the dip vector p represents the apparent 

dip of a reflector in the migrated domain, this reflector 
passing through the image point x . 

The dip vector p is expressed according to: 

p = {dx 3 /dxi, dx 3 /dx 2f 1) Equation (7a) 

2 0 where 5x 3 /9xi represents the apparent dip along the 

xi axis and dx 3 /dx 2 represents the apparent dip along the 
x 2 axis. 

0 and \|/ represent the inclination and azimuth 
respectively of the dipped reflector passing through the 

2 5 image point 3c at which the illumination is evaluated. 

Thus, the following expressions are obtained: 
cos9 = (l+Pi+p 2 )" 1/2 Equation (7b) 

tany = p 2 /pi Equation (7c) 

As we have seen above, starting from the Taylor 

30 series development of the diffraction travel time (known 
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10 X r = X + M . X p = X + 



Equation (7d) 



from the isochronic migration map) , the following 
equation is applicable: 

td(jfr) = t d (JC) + (V x t d ) T (X r -i)+^(i r -X ) T .A X(X t d . (5r 

3c ) + ... Equation (2) 

The reflection point x r is firstly estimated along 
the reflection plane for which the diffraction travel 
time is stationary (condition for specular reflection) . 

In order to consider positions only belonging to the 
dipped plane, a new vector x p is introduced such that: 

f\ 0 ^ 

0 1 
K~P\ -Pi. 

Note that in this case, as required, the two vectors 
of the matrix M are actually perpendicular to the dip 
vector p • 

Rewriting equation (2) using equation (7d) gives the 
15 following expression: 

t d (x p ) = t d (x) + (V x t d ) T M. x p + i .'x T p .M T .A X/X t d . Mx T p 

+... Equation (7e) 

namely T d (5c p ) = T d {x) + b T .x. p + ]^.x T p .Fx T p + ... 

Equation (7f) 



fx ^ 



20 In this equation (7f ) : 

b T = (V >: t d ) T M r where ■ b = M T .(V x t d ), and F = M T 
A x , x t d .M. 

Consequently, the position x rp for which the 
diffraction travel time-t d (x p ) is stationary is deduced: 
25 Vx P t d (5p) | xrp = 0 T => x^-F^.b =>■ x x =x+M.x rp 

Equation (7g) 
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Equations (7g) and (2) are combined to obtain an 
expression for the reflection travel time in equation 

(4) : t E = U(x z ) w U(x) - ^ fc'.F" 1 . b . 

As already mentioned, the illumination fold is 
5 incremented during operation 3 as a function of the 
difference between the diffraction travel time t d and the 
reflection travel time t r . 

Even if the reflection travel time t v (x r (p); s, r) 
is different from the diffraction travel time 
0 t d (x;s,r) f due to the limited band nature of seismic 
data, the migration of the reflection that took place at 
x can constructively contribute to the image of the 
reflector with a dip along x . 

The selection between illuminating migration 
5 contributions (constructive interference) and non- 
illuminating migration contributions (destructive 
interference) is made using the increment function 

i (t d , t r ; s , r ) . 

According to one embodiment of the invention, the 
0 increment function is expressed according to i(t d ,t r ; 
s,r) = s(t d (i;j,r) - t x (x r [p); s, r)) f where: 

- s(t) is a seismic wavelet which may be 
monochromatic, with limited band or infinite band and it 
can be chosen to be dependent on time and space, 
5 - s(t) is the derivative with respect to time of the 

seismic wavelet s(t) = d[s (t) ]/dt . 

A derivative with respect to time is usually applied 
to each seismic trace before the Kirchhoff migration (in 
time or in depth). Thus, particularly in order . to make 
0 the studied quantities uniform, the derivative s(t) of 
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the seismic wavelet with respect to time is chosen 
instead of using the seismic wavelets s(t) itself. 

According to another embodiment of the invention, 
particularly when the derivatives of seismic traces with 
respect to time are not determined, the increment 
function takes account of the seismic wavelet s(t) 
instead of the derivative s (t) of the wavelet with 
respect to time; the result is then i(t d ,t r ; s , r) = 
s(t 6 (x;s,r) - t r (x r (p ) ; s , f ) ) . 

A migration amplitude correction term w ( 3c ; s , r ) is 
conventionally calculated during the migration, and in 
particular takes account of a regular acquisition 
geometry . 

Therefore, this term w{x;s f r) may advantageously 
be used in the method according to the invention, and act 
as "preconditioning" to the correction made by 
determination of the previously described illumination 
fold. 

The increment function may then be weighted by this 
preconditioning term, such that the illumination fold 
I ( x , p ; s , r ) ) related to a (source s, receiver r) pair 
is incremented by: 

I(t d ,t r ; s , r ) . w ( x ; s , r ) . 

The evaluation of the illumination fold in the 
migrated domain according to the invention is optimum 
since it is applicable to any acquisition geometry. 

Based on an inverse approach, the illumination fold 
I(x,p) estimated by the invention agrees with the 
theoretical results established by Bleistein. This is not 
the case for other methods proposed to solve a similar 
problem. 
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As already mentioned earlier, the evaluation of the 
illumination fold proposed by the invention uses only 
basic quantities available for migration conventionally 
done by a Kirchhoff summation. 
5 Advantageously, the illumination fold according to 

the invention may be estimated while seismic traces are 
being migrated. 

This is not the case for some other proposed methods 
(Schneider et al . , Laurain and Vinje) for which the said 
10 evaluation is done post-migration. 

According to another aspect, the invention also 
proposes a method for correction of the amplitudes of 
seismic data recorded during 3D .seismic prospecting in 
order to compensate for the effect of non-uniform 
5 illumination of subsoil reflectors, characterized in that 
it comprises steps consisting of: 

- estimating the illumination fold I(x,p) by the 
method according to the first aspect of the invention 
described in detail above, 

0 - using the inverse I" 1 ^,^) of the said ratio as 

the weighing factor to be applied to each of the said 
seismic data amplitudes. 

According to yet another aspect, the invention also 
proposes a method for selection of one of a plurality of 
5 acquisition geometries as a function of the target of 3D 
seismic prospecting, characterized in that it comprises 
steps consisting of: 

- for each acquisition geometry considered, 
determining the illumination fold I ( x , p ) estimated 

0 according to the first aspect of the invention given in 
detail above, 
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- selecting the acquisition geometry providing the 
optimum illumination fold as a function of the target. 

Finally, it will be noted that the illumination fold 
I(x,p) estimated according to the first aspect of the 
invention described in detail above, may also be used 
after the migration to eliminate acquisition . patterns, 
and to deconvolute the migrated image. 
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22 
CLAIMS 

1. Method for estimating the seismic illumination 
fold \{x,p) in the migrated 3D domain at at least one 
image point x f for at least one dip of vector p 

wherein the illumination fold I(x,p;s,?) for each (source- ^ , 
5 receiver f ) pair in the seismic survey is estimated, by 

applying the following steps: 

- determination of the reflection travel time 

t r ( x r ( P ) ; s , r ) from the source s to the specular 

reflection point x r on the plane reflector ' passing 
10 through the image point x and perpendicular to the dip 

vector p and then returning to the reflector r ; 

starting from the diffraction travel time ta{x;s,r) 

from the source s to the said image point x and then 

returning to the reflector r ; 
.15 - incrementing the said illumination fold I(x,p;s,r) 

related to the said (source s r receiver r) pair as a 

function of the difference between the diffraction travel 

time -t d (x;«s 5 r) and the reflection travel time 

t r ( i r ( P ) ; I , r ) . 

20 2. Method according to claim 1, comprising the step 

of summating each of the said illumination folds I(x t p;s,?j 

related to a (source s , receiver r) pair so as to 

determine the total illumination fold I(x,p) = y i(x,p;ij; . 

. s,p 

3. Method according to one of the preceding claims, 
25 wherein, during the incrementing step, the illumination 
fold I(x i p,s i r) is incremented using an increment function 
i(t d ,t r ; s,r ) according to I(x.p) = I (x,p> + i(t d ,t r ; s,r ) , the 
said increment function taking account of the difference 
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between the diffraction travel time t d ( x ; s, r) and the 
reflection travel time t T {x z (p); s, r). 

4. Method according to claim 3, wherein the 
increment function i is a function of the seismic wavelet 
s(t). 

5. Method according to claim 4, wherein the 
increment function i is expressed as a function of the 
derivative of the seismic wavelet s(t) according to: 
i(t d/ t r ; s , r) = s(t d (x; s , r) - t r (x t (p);s,r) 

■ 6. Method according to claim 4, wherein the 
increment function i is expressed as a function of the 
derivative s (t) of the seismic wavelet s (t) with respect 
to time according to: 

i(t d , t r ; s , f ) = t d (x; s , r ) - t r ( 3r r ( p ) ; s , r ) 

7. Method according to any one of claims 3 to 6, in 
which an a priori correction w(x,s,r) of the illumination 
fold is taken into account by migration, comprising the 
step of incrementing the illumination fold I ( 3c , p ; s , r ) 
related to a (source s , receiver r) pair by 
i (t d/ t r ; s,r).vi(x;s,r) . 

8 . Method according to any one of the preceding 
claims, wherein the determination step includes the 
second order Taylor series development of the diffraction 
travel time t d (x;^,F) around the image point x : 

t d < X ; S , T )=t d ( X ; S , f ) + ( V «t d { X ; S , F ) ) T . ( X r -X )+^ ( X r ~ 
X) T .A, :/X t d (X; s , r ) . (X r -X) 

9. Method according to claim 8, wherein the specular 

reflection point x r ( P ) is determined along the length of 
the said reflector such that the diffraction travel time 

at the said specular reflection point x * (P ) is 
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stationary, according to the equation: 

p t A( V x t d { x; s , r )+A,, ;: t d ( x ; s , r ) . ( x r ( P ) -x ) )-o . 

10. Method according to any one of claims 8 or 9, 
wherein the specular reflection point x r and the 

5 reflection travel time t r ( x r (p);s r r ) (are determined 
according to the following expressions: 

x r (p) = x-M.F -1 . b 

t r {Xr(p);s r r) = t d ( x ; s , r ) . b r . F^ 1 . b 

10 

where: 

- M is a (3x2) matrix described by two vectors 
extending along the length of the reflector, and 
therefore perpendicular to the dip vector p ; 

15 - b is a (2 x 1) vector of first order derivatives 

of the diffraction travel time along the reflection 

plane: b = M T . ( V x t d ) ; 

• - F is a (2 x 2) matrix of second order derivatives 
of the diffraction travel time along the reflection 
20 plane: F = M T . (A x , x .t d ) .M. 

11. Method according to any one of the preceding 
claims, wherein the determination step uses isochronic 
migration maps t d ( x;?.;" ) specified for each (source s f 
receiver r ) pair involved in the migration at each image 

25 point x in the migrated 3D domain. 

12. Method according to any one of the preceding 
• claims, wherein the seismic illumination fold I(x,p) in 

the migrated 3D domain is estimated during the Kirchhoff 
summation migration of seismic data recorded during the 
30 3D seismic prospecting. 
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13. Method for correction of seismic data amplitudes 
recorded during 3D seismic prospecting in order to 
compensate for the effect of non-uniform illumination of 
sub-soil reflectors, comprising the steps of: 
5 - estimating the illumination fold I(x 9 p) using the 

method according to any one of claims 1 to 12, 

- using the- inverse l' x (x,p) of the said ratio as a 
weighting factor to be applied to each of the said 
seismic data amplitudes. 

10 14. Method for selection of an acquisition geometry 

among a plurality of acquisition geometries as a function 
of the target of 3D seismic prospecting, comprising the 
steps of: 

- determining the illumination fold I{x,p) by the 
15 method according to any one of claims 1 to 12, for each 

of the acquisition geometries considered, 

- selecting the acquisition geometry providing the 
optimum illumination fold as a function of the target. 
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For each (source s, receiver?*) pair involved in the migration. 



For each image point x in the migrated domain. 



For each dip component p selected for the illumination 
calculation I (xfp,!;, T). 



Determine the reflection travel time tr (xr (p ) ; s, r) at the 
specular reflection point x^on the plane reflector passing 
through the image point x*and perpendicular to the dip vector ~p. 



Increment the illumination fold I (>Cp,%T) relative to the 
(source receiver/) pair as a function of the difference between: 

the diffraction travel time (t d (x*; %T) and 

the reflection travel time t r (x? (p ) ; s| T). 
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